### Make sure all the packages are installed
if (!requireNamespace("Seurat", quietly = TRUE))
install.packages("Seurat")
if (!requireNamespace("tidyverse", quietly = TRUE))
install.packages("tidyverse")
if (!requireNamespace("colorBlindness", quietly = TRUE))
install.packages("colorBlindness")
if (!requireNamespace("RColorBrewer", quietly = TRUE))
install.packages("RColorBrewer")
if (!requireNamespace("DT", quietly = TRUE))
install.packages("DT")
if (!requireNamespace("ComplexHeatmap", quietly = TRUE))
BiocManager::install("ComplexHeatmap")
if (!requireNamespace("openxlsx", quietly = TRUE))
BiocManager::install("openxlsx")
if (!requireNamespace("presto", quietly = TRUE))
devtools::install_github("immunogenomics/presto")
if (!requireNamespace("SeuratData", quietly = TRUE))
devtools::install_github('satijalab/seurat-data')
### Load all the necessary libraries
library(Seurat)
library(tidyverse)
library(colorBlindness)
library(RColorBrewer)
library(DT)
library(ComplexHeatmap)
library(SeuratData)
library(openxlsx)7 - Differential Gene Expression & Level 1 Annotation
Introduction
In this notebook we are going to pick up from the clustering notebook and focus on how to annotate a single-cell dataset. We’re briefly going to touch on how to carry out differential expression between clusters but for a more in depth explanation visit our notebook from a previous workshop. Annotating single-cell data is the most laborious part of the analysis. It requires in depth knowledge of the cell types making up your biological system and multiple rounds of iterative clustering and subsetting to get to a fine grained annotation. Luckily, there are some tools, like label transfer, that can be used to help speed up this process. These tools rely on using a relevant and previously annotated reference to help automate this annotation process. These automated methods work great for coarse cell type labels but don’t perform quite as well when annotating fine grained cell states that may be not be found in our reference. Therefore, an approach combining both methodologies is usually used. A good resource going over this process is the paper by Luecken, M, et al..
Key Takeaways
To annotate our clusters we need to determine which genes are differentially expressed in each one.
Differentially expressed genes depend on which cell types we are comparing. The same cell type will have different differentially expressed genes if we change the other cell types in the dataset.
We can quantify these differentially expressed genes using effect size and discriminatory power metrics such as log2FC and AUC.
P values obtained from carrying out DGE analysis between clusters are inflated and should not be used.
Annotation is a laborious process that uses automated and manual approache in which we use automated methods to coarsley label cells and manual annotation to identify unique cell states.
Automated annotation requires relevant reference data that has been previously annotated. We will only be able to identify cell types that were anotated in the reference.
Manual annotation is based on literature knowledge and digging through previous annotation efforts or field-relevant papers describing cell types of interest using other methods - bulk RNAseq, FACS…
Libraries
Load data
We’re going to be working with a dataset from the paper - Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19 Download data from cellxgene portal.
# Download the data in data/ directory
# download.file(
# url = "https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds",
# destfile = "../data/workshop-data.rds",
# method = "wget",
# extra = "-r -p --random-wait")
# We can also use the CLI with the wget command below
# wget https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds
se <- readRDS("../data/d8e35450-de43-451a-9979-276eac688bce.rds")Generate a color palette for our cell types
# https://www.datanovia.com/en/blog/easy-way-to-expand-color-palettes-in-r/
# nb.cols <- length(unique(se$Celltype))
# mycolors <- colorRampPalette(paletteMartin)(nb.cols)
pal <- paletteMartin
names(pal) <- sort(unique(se$Celltype))Analysis
Convert ENSEMBL IDs to Gene Symbols
Right away we can see how ensembl ids are used in the rownames. Let’s transform them into their matched symbols to make them human-readable:
head(rownames(se))[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" "ENSG00000000938"Convert to gene symbols
gene_df <- readr::read_csv(file = "../data/cov_flu_gene_names_table.csv")
symbol_id <- data.frame(index = rownames(se)) %>%
left_join(gene_df, by = "index") %>%
pull(feature_name)
# re-create seurat object
mtx <- se@assays$RNA$data
rownames(mtx) <- symbol_id
se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
rm(mtx); gc() used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5392750 288.1 10451644 558.2 NA 10420433 556.6
Vcells 183647299 1401.2 674232762 5144.0 98304 821666950 6268.9Quick processing
se <- se %>%
NormalizeData(verbose = FALSE) %>%
FindVariableFeatures(
method = "vst",
nfeatures = 3000,
verbose = FALSE) %>%
ScaleData(verbose = FALSE, features = VariableFeatures(.)) %>%
RunPCA(verbose = FALSE)
ElbowPlot(se, ndims = 50)
Let’s run FindNeighbors and FindClusters to label our data:
se <- FindNeighbors(se, reduction = "pca") %>%
FindClusters(resolution = c(0.01, 0.05, 0.1, 0.25))Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9956
Number of communities: 4
Elapsed time: 10 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9829
Number of communities: 7
Elapsed time: 11 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9721
Number of communities: 9
Elapsed time: 10 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9532
Number of communities: 15
Elapsed time: 12 secondsLet’s compute the UMAP
se <- se %>%
RunUMAP(dims = 1:30, verbose = FALSE)We can now take a look on the UMAP space how the clusters look
se$sample_id <- se$`Sample ID`
DimPlot(
se,
group.by = c(
"RNA_snn_res.0.01", "RNA_snn_res.0.05",
"RNA_snn_res.0.1", "RNA_snn_res.0.25"),
label = TRUE)
dim_plt <- DimPlot(
se,
group.by = c("RNA_snn_res.0.05"),
label = TRUE)And the original cell type labels + the sample IDs
DimPlot(
se,
group.by = c("Celltype", "sample_id"),
label = FALSE)
For the purpose of this tutorial we’re going to go forward with resolution 0.05!
DGE Wilcoxon
The different implementations Seurat incorporates provides in FindAllMarkers compare the gene expression between 2 groups of cells. This one vs all strategy is very quick and returns the avg_log2FC. This avg_log2FC is computed as detailed here & here. Since we’re working with normalized data the log2FC can be directly computed by subtracting the average expression between both groups - log_{2}(\frac{exp1}{exp2})=log_{2}(Avg\_exp1)-log_{2}(Avg\_exp2)
Idents(se) <- "RNA_snn_res.0.05"
mgs <- FindAllMarkers(
se,
test.use = "wilcox",
slot = "data",
only.pos = TRUE,
logfc.threshold = 0.5,
min.pct = 0.25)Look at the results in a dynamic table:
DT::datatable(mgs, filter = "top")Look at the results in a heatmap
top10 <- mgs %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1) %>%
slice_head(n = 10) %>%
ungroup()
DoHeatmap(se, features = top10$gene) + NoLegend()
Save marker genes to a spreadsheet
mgs_ls <- lapply(unique(mgs$cluster), function(i) {
mgs %>% dplyr::filter(cluster == i)
})
# Set names which will be the sheet name
names(mgs_ls) <- unique(mgs$cluster)
# Save marker genes to spreadsheet
openxlsx::write.xlsx(mgs_ls, file = "../data/lvl1_mgs.xlsx")P-value interpretation
More details can be obtained in OSCA.
P values obtained from DGE analysis are inflated and, therefore invalid in their interpretation. We can’t use p-values to reject the Null Hypothesis since we are carrying out data snooping. This means that we are dividing the clusters based on their gene expression, and then computing p-values for the genes that are differentially expressed, even though we know these clusters have different gene expression patterns since we clustered the data based on them being different.
A way to show this is by looking at how skewed the distributions of the p-values obtained is:
# Compute the p-values without he thresholds
mgs2 <- FindAllMarkers(
se,
test.use = "wilcox",
only.pos = TRUE,
logfc.threshold = 0,
min.pct = 0,
return.thresh = 1,
max.cells.per.ident = 100 # use 100 cells per cell type for speed
)
ggplot(mgs2, aes(x = p_val, fill = cluster, color = cluster)) +
# geom_histogram(alpha = 0.3, position = "identity") +
geom_density(alpha = 0.3) +
theme_minimal()
ggplot(mgs2, aes(x = p_val, fill = cluster, color = cluster)) +
geom_histogram(alpha = 0.3, position = "identity") +
facet_wrap(~cluster, scales = "free") +
theme_minimal()
Annotation
There are two main ways to annotate your single cell dataset which are usually used together to aid in the process:
- Automatic cell type annotation: typicall requires a reference dataset representative of the samples we are trying to annotate. Some commonly used tools in R are
SingleRwhich has very good documentation andSeurat’s reference mapping either through their Azimuth server or following their reference mapping vignette which performs remarkably well.
Azimuth
We can use Azimuth online app to manually use a reference dataset to annotate our dataset.
To do so we need to save the counts matrix as an RDS file.
se@assays$RNA$counts[1:15, 1:15]15 x 15 sparse Matrix of class "dgCMatrix"
TSPAN6 . . . . . . . . . . . . . . .
TNMD . . . . . . . . . . . . . . .
DPM1 . . . . . . . 1 . 1 . . 1 . .
SCYL3 . . . . . . . . . . . . . . .
C1orf112 . . . . . 1 . . . . . . . . .
FGR . . 1 2 . . 3 . 1 1 1 . 2 . .
CFH . . . . . . . . . . . . . . .
FUCA2 . . . . . . 1 . . . . . . . .
GCLC . . . . 3 . 2 . . 1 . . . . .
NFYA . . . . . . . . . . . . . . .
STPG1 . . . . . . . . . . . . . . .
NIPAL3 . 1 . . . . . . . . . . . . .
LAS1L . . . . . . . . . . . . . . .
ENPP4 . . . . . . . 1 . . . 1 . . .
SEMA3F . . . . . . . . . . . . . . .saveRDS(
object = se@assays$RNA$counts,
file = "../data/counts.rds")Reference Mapping
We are going to download a reference human PBMC dataset using the SeuratData package. You can use whichever reference suits your dataset best.
SeuratData::AvailableData() Dataset Version Summary species system ncells tech seurat default.dataset disk.datasets other.datasets notes Installed InstalledVersion
adiposeref.SeuratData adiposeref 1.0.0 Azimuth Reference: adipose human adipose 160075 scRNA-seq and sNuc-seq <NA> <NA> <NA> <NA> <NA> FALSE <NA>
bmcite.SeuratData bmcite 0.3.0 30k Bone Marrow Cells human bone marrow 30672 <NA> 3.2.2 <NA> <NA> <NA> <NA> FALSE <NA>
bonemarrowref.SeuratData bonemarrowref 1.0.0 Azimuth Reference: bonemarrow human bonemarrow 297627 10x v2 <NA> <NA> <NA> <NA> <NA> FALSE <NA>
cbmc.SeuratData cbmc 3.1.4 scRNAseq and 13-antibody sequencing of CBMCs human CBMC (cord blood) 8617 CITE-seq 3.1.4 raw processed <NA> <NA> FALSE <NA>
celegans.embryo.SeuratData celegans.embryo 0.1.0 6k C. elegans embryos from Packer and Zhu et al (2019) C. elegans embryo 6188 <NA> <NA> raw <NA> <NA> <NA> FALSE <NA>
fetusref.SeuratData fetusref 1.0.0 Azimuth Reference: fetus human fetus 377456 <NA> <NA> <NA> <NA> <NA> <NA> FALSE <NA>
hcabm40k.SeuratData hcabm40k 3.0.0 40,000 Cells From the Human Cell Atlas ICA Bone Marrow Dataset human bone marrow 40000 10x v2 <NA> raw <NA> <NA> <NA> FALSE <NA>
heartref.SeuratData heartref 1.0.0 Azimuth Reference: heart human heart 656509 scRNA-seq and sNuc-seq <NA> <NA> <NA> <NA> <NA> FALSE <NA>
humancortexref.SeuratData humancortexref 1.0.0 Azimuth Reference: humancortex human motor cortex 76533 <NA> <NA> <NA> <NA> <NA> <NA> FALSE <NA>
ifnb.SeuratData ifnb 3.1.0 IFNB-Stimulated and Control PBMCs human PBMC 13999 10x v1 <NA> raw <NA> processed <NA> FALSE <NA>
kidneyref.SeuratData kidneyref 1.0.2 Azimuth Reference: kidney human kidney 64693 snRNA-seq <NA> <NA> <NA> <NA> <NA> FALSE <NA>
lungref.SeuratData lungref 2.0.0 Azimuth Reference: lung human lung 584884 <NA> <NA> <NA> <NA> <NA> <NA> FALSE <NA>
mousecortexref.SeuratData mousecortexref 1.0.0 Azimuth Reference: mousecortex mouse motor cortex 159738 10x v3 <NA> <NA> <NA> <NA> <NA> FALSE <NA>
panc8.SeuratData panc8 3.0.2 Eight Pancreas Datasets Across Five Technologies human Pancreatic Islets 14892 SMARTSeq2, Fluidigm C1, CelSeq, CelSeq2, inDrops <NA> raw <NA> <NA> <NA> TRUE 3.0.2
pancreasref.SeuratData pancreasref 1.0.0 Azimuth Reference: pancreas human pancreas 35289 <NA> <NA> <NA> <NA> <NA> <NA> FALSE <NA>
pbmc3k.SeuratData pbmc3k 3.1.4 3k PBMCs from 10X Genomics human PBMC 2700 10x v1 3.1.4 raw <NA> pbmc3k.final <NA> TRUE 3.1.4
pbmcMultiome.SeuratData pbmcMultiome 0.1.4 10X Genomics PBMC Multiome Dataset human pbmc 11909 <NA> 4.1.1 NA <NA> pbmc.rna, pbmc.atac One sample with two modalities, RNA and ATAC FALSE <NA>
pbmcref.SeuratData pbmcref 1.0.0 Azimuth Reference: pbmc human PBMC 2700 10x v1 <NA> <NA> <NA> <NA> <NA> TRUE 1.0.0
pbmcsca.SeuratData pbmcsca 3.0.0 Broad Institute PBMC Systematic Comparative Analysis human PBMC 31021 10x v2, 10x v3, SMARTSeq2, Seq-Well, inDrops, Drop-seq, CelSeq2 <NA> raw <NA> <NA> HCA benchmark TRUE 3.0.0
ssHippo.SeuratData ssHippo 3.1.4 Slide-seq v2 dataset of mouse hippocampus mouse hippocampus NA slideseq v2 <NA> raw <NA> <NA> <NA> FALSE <NA>
stxBrain.SeuratData stxBrain 0.1.1 10X Genomics Visium Mouse Brain Dataset mouse brain 12167 visium <NA> NA <NA> posterior1, posterior2, anterior1, anterior2 One sample split across four datasets as paired anterior/posterior slices FALSE <NA>
stxKidney.SeuratData stxKidney 0.1.0 10X Genomics Visium Mouse Kidney Dataset mouse kidney 1438 visium <NA> raw <NA> <NA> <NA> FALSE <NA>
thp1.eccite.SeuratData thp1.eccite 3.1.5 ECCITE-seq THP-1 human <NA> NA <NA> <NA> <NA> <NA> <NA> <NA> FALSE <NA>
tonsilref.SeuratData tonsilref 2.0.0 Azimuth Reference: tonsil human tonsil 377963 scRNA-seq <NA> <NA> <NA> <NA> <NA> FALSE <NA>SeuratData::InstallData("pbmcsca", force.reinstall = TRUE)
SeuratData::InstalledData() Dataset Version Summary species system ncells tech seurat default.dataset disk.datasets other.datasets notes Installed InstalledVersion
panc8.SeuratData panc8 3.0.2 Eight Pancreas Datasets Across Five Technologies human Pancreatic Islets 14892 SMARTSeq2, Fluidigm C1, CelSeq, CelSeq2, inDrops <NA> raw <NA> <NA> <NA> TRUE 3.0.2
pbmc3k.SeuratData pbmc3k 3.1.4 3k PBMCs from 10X Genomics human PBMC 2700 10x v1 3.1.4 raw <NA> pbmc3k.final <NA> TRUE 3.1.4
pbmcref.SeuratData pbmcref 1.0.0 Azimuth Reference: pbmc human PBMC 2700 10x v1 <NA> <NA> <NA> <NA> <NA> TRUE 1.0.0
pbmcsca.SeuratData pbmcsca 3.0.0 Broad Institute PBMC Systematic Comparative Analysis human PBMC 31021 10x v2, 10x v3, SMARTSeq2, Seq-Well, inDrops, Drop-seq, CelSeq2 <NA> raw <NA> <NA> HCA benchmark TRUE 3.0.0ref <- LoadData("pbmcsca")Look at the data
table(ref$CellType)
B cell CD14+ monocyte CD16+ monocyte CD4+ T cell Cytotoxic T cell Dendritic cell Megakaryocyte Natural killer cell Plasmacytoid dendritic cell Unassigned
5020 5550 804 7391 9071 433 977 1565 164 46 Preprocess data - Normalize, identify HVG and look at the Elbow plot. We start by normalizing the data just to make sure both datasets are processed the same way. We then use the union of the HVG to make sure we capture the biological variability of both datasets. Next we compute the UMAP, it is necessary to specify return.model=TRUE! Lastly, we look at the Elbow plot to assess how many dimensions to use in downstream analysis.
# Renormalize data and find HVG to make sure
se <- NormalizeData(se) %>% FindVariableFeatures(method = "vst", nfeatures = 3000)
ref <- NormalizeData(ref) %>%
FindVariableFeatures(method = "vst", nfeatures = 3000) %>%
ScaleData() %>%
RunPCA()
ElbowPlot(ref, ndims = 30)
ref <- ref %>% Seurat::RunUMAP(reduction = "pca", dims = 1:30, return.model = TRUE)
# Use the union of HVG
hvg <- union(VariableFeatures(se), VariableFeatures(ref))
length(hvg)[1] 5464Reference mapping - there is very good explanation in the documentation of ?FindTransferAnchors. Basically what we are doing here is embedding the cells into the same latent space. Anchors - defined as pairs of cells contained within each other’s k-neighborhood - are identified. Lastly, with TransferData the shared nearest neighbors overlap between the anchors and the cells is used to obtain the that cell’s predicted label.
# Find anchors in reference dataers
anchors <- FindTransferAnchors(
reference = ref,
query = se,
dims = 1:30,
reference.reduction = "pca",
features = hvg,
normalization.method = "LogNormalize")
# Extract the predictions
predictions <- TransferData(
anchorset = anchors,
refdata = ref$CellType,
dims = 1:30)
head(predictions) predicted.id prediction.score.Cytotoxic.T.cell prediction.score.CD4..T.cell prediction.score.CD14..monocyte prediction.score.B.cell prediction.score.Megakaryocyte prediction.score.CD16..monocyte prediction.score.Natural.killer.cell prediction.score.Dendritic.cell prediction.score.Plasmacytoid.dendritic.cell prediction.score.Unassigned prediction.score.max
AAACCCAAGAATGTTG-12 Natural killer cell 0.21689575 0 0 0 0 0 0.78310425 0 0 0 0.7831042
AAACCCAAGCATTGTC-19 CD4+ T cell 0.00000000 1 0 0 0 0 0.00000000 0 0 0 1.0000000
AAACCCAAGCTACGTT-6 Cytotoxic T cell 0.98981523 0 0 0 0 0 0.01018477 0 0 0 0.9898152
AAACCCAAGGCCGCTT-3 CD14+ monocyte 0.00000000 0 1 0 0 0 0.00000000 0 0 0 1.0000000
AAACCCAAGGCCTGCT-12 Natural killer cell 0.08903623 0 0 0 0 0 0.91096377 0 0 0 0.9109638
AAACCCAAGGGCAATC-1 B cell 0.00000000 0 0 1 0 0 0.00000000 0 0 0 1.0000000Let’s add these predictions to our seurat object
colnames(predictions) [1] "predicted.id" "prediction.score.Cytotoxic.T.cell" "prediction.score.CD4..T.cell" "prediction.score.CD14..monocyte" "prediction.score.B.cell" "prediction.score.Megakaryocyte" "prediction.score.CD16..monocyte" "prediction.score.Natural.killer.cell" "prediction.score.Dendritic.cell" "prediction.score.Plasmacytoid.dendritic.cell" "prediction.score.Unassigned" "prediction.score.max" se <- AddMetaData(
se,
metadata = predictions[, c("predicted.id", "prediction.score.max", "prediction.score.Unassigned")])Look at the label transfer
DimPlot(se, group.by = c("Celltype", "predicted.id"))
Since we also got a confidence value for each cell we can visualize it
FeaturePlot(se, features = c("prediction.score.max", "prediction.score.Unassigned"))
Manual Annotation
Cluster 0 & 4
Let’s look at genes that are differentially expressed
FeaturePlot(
se,
features = c("CD3D", "CD3D", "TRAC", "TRBC2", "CD8B", "CD4")) +
dim_plt
VlnPlot(
se,
features = c("CD3D", "CD3D", "TRAC", "TRBC2", "CD8B", "CD4"),
group.by = "RNA_snn_res.0.05") +
dim_plt
Clusters 0 & 4 seem to have a lot of expression of T cell related genes so at this level 1 we are going to label them as T cells.
Cluster 1
Let’s look at genes that are differentially expressed
FeaturePlot(
se,
features = c("CD14", "S100A8", "VCAN", "LYZ", "MS4A6A")) +
dim_plt
VlnPlot(
se,
features = c("CD14", "S100A8", "VCAN", "LYZ", "MS4A6A"),
group.by = "RNA_snn_res.0.05") +
dim_plt
Cluster 1 is expressing a lot of monocyte-like genes, at this level 1 we are going to label them as monocytes.
Cluster 2
Let’s look at genes that are differentially expressed
FeaturePlot(
se,
features = c("MS4A1", "CD79A", "CD79B", "IGHD", "IGHM")) +
dim_plt